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Energy landscape mappings are performed for two different molecular systems under mechanical 
loads. Barrier heights are observed to scale as AU ~ S^^^, where 5 is a residual load. Catastrophe 
theory predicts that this scaling should arise for vanishing S, however, this region is irrelevant in 
physical processes at finite temperature because thermal fluctuations cause the system to cross over 
the barrier before reaching the small S regime. Surprisingly, we find that the AU ~ S^^^ scaling is 
valid far beyond the vanishing S regime described by catastrophe theory. This scaling will therefore 
be relevant at finite temperatures, and can be the basis for corrections to standard rate theoretic 
approaches. 



In a broad range of condensed matter systems, one is 
interested in the question of how some material responds 
to an external mechanical load. External loads cause 
liquids to flow, in Newtonian or various types of non- 
Newtonian flows. Glassy materials, composed of poly- 
mers, metals, or ceramics, can deform under mechanical 
loads, and the nature of the response to loads often dic- 
tates the choice of material in various industrial applica- 
tions. In biological systems, the response of proteins to 
external loads governs aspects of cell adhesion and muscle 
function. 

The nature of all of these responses depends on both 
the temperature and loading rate. As described by 
Eyring mechanical loading lowers energy barriers, 
thus facilitating progress over the barrier by random ther- 
mal fluctuations. The Eyring model approximates the 
loading dependence of the barrier height as linear. The 
Eyring model, with this linear barrier height dependence 
on load, has been used over a large fraction of the last 
century to describe the response of a wide range of sys- 
tems and underlies modern approaches to biophysical 
rupture processes 0, 0, 13 j sheared glasses 0, , etc. 

The linear dependence will always correctly describe 
small changes in the barrier height, since it is simply the 
first term in the Taylor expansion of the barrier height 
as a function of load. It is thus appropriate when the 
barrier height changes only slightly before the system es- 
capes the local energy minimum. This situation occurs 
at higher temperatures; for example, Newtonian flow is 
obtained in the Eyring model in the limit where the sys- 
tem experiences only small changes in the barrier height 
before thermally escaping the energy minimum. 

As the temperature decreases, larger changes in the 
barrier height occur before the system escapes the energy 
minimum (giving rise to, for example, non-Newtonian 
flow). In this regime, the linear dependence is not neces- 
sarily appropriate, and can lead to inaccurate modeling. 
For example, Li and Makarov have shown that there 
is a nonlinear barrier height dependence in stretched pro- 
teins, and that the assumption of a linear dependence in 
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FIG. 1: Left: Energy from equation |2| for various S. 5 — 0: 
solid, S = -1/8: dotted, 6 = -1/4: dashed. Right: Bifur- 
cation diagram indicating the locations of the extrema as a 
function S. Minima: solid, Barrier: dotted. At 5 = the min- 
imum of interest collides with a barrier. At 5 — Smin ~ —1/4, 
the barrier collides with a distant minimum and ceases to ex- 
ist; it makes no sense to discuss quantities such as AU for 

the analysis of experimental results leads to inaccurate 
conclusions. 

The present investigation addresses this load depen- 
dence of the barrier height. The analysis is based on the 
energy landscape formalism which considers dynam- 
ics to be the sum of vibrational- like motion within energy 
minima and transitions between energy minima. Barriers 
are associated with saddle points that connect adjacent 
energy minima. 

In molecular systems, the energy is a smooth function 
of the internal degrees of freedom plus a control parame- 
ter {e.g. the controlled stress or strain). As the con- 
trol parameter is varied, any minimum in the landscape 
will flatten out in some direction as the minimum collides 
with a first order saddle point collides [ifil Ull (c./. 
figure ^1 . This type of externally induced topolo gica l 
change in a function is known as a fold catastrophe [iSj. 

It has long been appreciated that these fold catastro- 
phes induce universal scalings of particular features of 
the surface in the limit where the minimum and saddle 
point are infinitesimally close together In this limit, 
the function has the lowest order Taylor expansion 

U = -Ax^ - Bx5 (1) 
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because the first order term (in the internal degrees of 
freedom) is zero at a minimum or saddle point, and the 
second order term (projected along the direction which 
connects the minimum and barrier) is zero at the point 
where the minimum and saddle point collide. In equa- 
tion^ X is the projection of the system's coordinates onto 
the zero curvature direction, S is the control parameter's 
distance away from the singularity, and A and B are pos- 
itive constants. To obtain the scaling laws, we note that 
the minimum and saddle point, X- and x+, are the points 
where the energy derivative vanishes, the curvature along 
the reaction coordinate at the minimum and saddle point 
correspond to the second partial derivatives of the func- 
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The fold ratio, rj = 6A[//(2A_a;'i), is unity when all 
three of these scaling relations are valid. UJ These ar- 
guments appeal only to (5's role as a control parameter 
and are equally valid when S represents, e.g.^ an imposed 
strain or stress. Recently, the consequences of these scal- 
ings have been discussed quantitatively in the context 
of phenomenological models [ia, llfil llTj . Fold ratios in 
incipient catastrophes have been measured in molecular 
simulations ,14] , but the individual scaling relations have 
not previously been addressed in externally driven molec- 
ular level simulations. 

These scaling relations must be obtained in the limit 
(5 — » 0, but at finite temperature, thermal fluctuations 
cause the barrier to be crossed before this vanishing 5 
regime is being reached. Little attention has been paid to 
the accuracy of these scaling relations for the physically 
meaningful finite S regime. While the fold ratio has been 
shown to deviate significantly from unity at finite S, 
the accuracy of the scaling relations for the individual 
quantities has not previously been addressed. 

Clues to how the scaling breaks down at finite \5\ are 
obtained from a l-flD energy function, U{x,6), based 
on arguments similar to those given in .11.] . The only 
requirement we make of U is that the coupling to the 
control parameter is bi-linear over the region of inter- 
est: U{x,S) — Uo{x) — BqxS. Demanding that (dU/dx) 
remain zero as we change 6, requires that 
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where xq is a stationary point. The energy of a stationary 
point then changes according to 
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where the second equality follows from mechanical equi- 
librium. The fold scalings obey the above relations be- 
tween the energy, position and curvature of a stationary 
point, but these latter relations are more general because 
they are based on much weaker assumptions about the 
form of the energy function than the cubic form used to 
obtain the fold scalings. Since the energy barrier is ob- 
tained after two integrations of the inverse curvature, we 
anticipate that as the load is backed away from the catas- 
trophe deviations from the scaling relations should occur 
first for the curvature then for the position and finally 
for the energy; i.e., the barrier height scaling should be 
more accurate at finite S than the other scaling relations. 

We first test this hypothesis regarding the relative ac- 
curacy of the barrier height scaling on a simple, analyt- 
ically solvable, model. This simple model includes the 
next order term in the Taylor expansion for the energy, 
giving a 4-th order polynomial. 



U = -x-' 



xS 
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The landscape for this energy function for various 6 is 
shown in figure . As S goes to zero, the (left hand) 

minimum and the energy barrier join together at a; = 0. 
As S is backed away from zero, the minimum and bar- 
rier move apart. As S is backed further away from zero, 
the barrier eventually collides with some other minimum 
at (5 = —1/4, rendering the quantities AU, X+, x+ un- 
defined. The results for this simple model are consistent 
with our expectation that the scaling of the barrier height 
will be more accurate than the scaling of the curvature: 
at the largest possible values of S, the barrier height scal- 
ing remains accurate to within 10%, while the scaling for 
A_ is in error by 100%, and the scaling prediction for A+ 
is infinite in error (A+ vanishes altogether). Analogous 
results are obtained with a negative fourth order term, in 
which case, the maximal \S\ occurs when the minimum 
collides with the other barrier (in this case A4. is at 50% 
of the value from the scaling relation, and A_ vanishes 
altogether). Similar results are obtained for other simple, 
analytically solvable, models. 

We have tested these ideas in simulations of realistic 
atomistic models by tracking local minima and saddle 
points as a control parameter is varied. To check the 
generality of our arguments, we investigate two very dif- 
ferent systems (a model glass, and a model protein), 
and consider both strain and stress as control param- 
eters. The model glass simulations use the Stillinger- 
Weber 80:20 mixture 0, and include 500 particles in 
an orthorhombic simulation cell with periodic boundary 
conditions (the Stillinger- Weber functional form, which 
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FIG. 2: From top to bottomiAC/, Ax, A- and A+ (plus sym- 
bols), and |r/ — 1| upon approach to a typical catastrophe. 
The solid (red) lines are the theoretical scaling predictions 
(equations, l^l 121 andEJ with prefactors determined by a fit to 
the smallest two decades of 7c — 7 shown. 7c is determined 
via optimization of At/ to the (7c — 7)'^'^^ form. 



is similar to the Lennard-Jones potential, is used because 
the continuous derivatives at a finite cutoff are necessary 
to analyze the energy landscape at the required preci- 
sion). The model protein is the Thirumalai model |l9l| . 
which consists of 46 sites that interact with bond stretch- 
ing, angle bending, torsion, and nonbonded (Lennard- 
Jones) interactions. A minimum is tracked by repeatedly 
minimizing the energy as the control parameter is varied 
in small increments. Similarly, a saddle point is tracked 
by repeatedly minimizing the sum of the squares of the 
forces as the control parameter is varied in small incre- 
ments (the saddle point is found initially by searching 
halfway along the vector that connects two minima). All 
numerical minimizations are performed using a variable 
metric algorithm, and the eigenvalues and eigenvectors of 
the system are computed using a standard QL reduction 
algorithm ^20]. 

Results for the glass, with shear strain as the control 




FIG. 3: Top: At/ and Bottom: A_ for the protein model in 
a length controlled mode as functions oi Lc — L where L is 
end to end length and Lc is the value of length at which the 
minimum and barrier collide. 

parameter, are shown in figure |21 Since the system is 
multidimensional, it will have many eigenvalues at the 
minimum, and we take A_ to be the smallest of these. 
The magnitude of the single negative eigenvalue at the 
barrier is denoted as A+ . We choose to use A_ to compute 
the fold ratio. 

In the small S limit all of the scaling relations (equa- 
tionsEllSl andEl are accurate, and the fold ratio is unity. 
This behavior is of course expected, because the catastro- 
phe is of the fold type. In terms of the large S behavior, 
the results are fully consistent with the ideas based on the 
arguments above, namely: (i) the barrier eigenvalue goes 
to zero at about 7c — 7 '^-^ .007 in a collision with some 
other minimum, (ii) the accuracy of the scaling relation 
for the curvature ( |2Jl is quite poor in comparison with 
the accuracy of the scaling relation for the barrier height 
(^, with the barrier height scaling being a reasonable 
approximation over the entire interval up to 7c— 7 ~ .007. 

For the protein model, the end to end distance of 
the protein, L, can be taken as the control parameter. 
As shown in figure |21 both increasing and decreasing L 
causes the lowest curvature at the minimum to go to 
zero, indicating the onset of catastrophes. As expected, 
the scaling relations are accurate at small values of S, 
where S, is Lc — L and Lc is the length at which the min- 
imum and barrier collide. At large S, the scaling relation 
for the curvature becomes poor while the scaling relation 
for the barrier height is reasonably accurate over the en- 
tire range of S. In contrast to the Lennard-Jones case, it 
is the minimum which disappears in a collision with an 
extraneous barrier with A_ going to zero at the maximal 
Lc — L, but, in both cases, the barrier scalings are found 
to be superior to the scalings of either of the curvatures. 

To use force as a control parameter in the protein 
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FIG. 4: Top: AI7 and Bottom: A_ for the protein model, as 
functions oi Fc — F where F is the external force and Fc is 
the value of force at which the minimum and barrier collide. 

model, L is reinstated as a bona fide degree of free- 
dom, and an external coupling is introduced so that 
?7tot = C^int ~ FL, where Uint is the usual internal energy 
of the protein. As the arguments for the fold scaling re- 
lations appeal only to the smoothness of the energy func- 
tion and not the particular mode of loading, we expect 
analogous results when force is the control parameter. 
AU and A_ are plotted in figure ^ Again, all scahng 
relations are accurate at small 6, but at large 6 the scal- 
ing relation for the AU is much more accurate than the 
scaling relation for A_. 

In summary, barrier heights in molecular systems are 
found to follow, to fairly high accuracy, the scaling re- 
lation AU ~ 5'^/^. While this scaling relation has been 
known to be rigorously valid in the (5 — > limit, this van- 
ishing S regime is not not physically significant at finite 
temperature, because thermal fluctuations cause the sys- 
tem to cross the barrier before the low delta regime is 
reached. However, our investigation shows that the scal- 
ing relation is appropriate outside of the low S limit - 
even when the scalings fail dramatically for the the cur- 
vatures or the fold ratio. The barrier height scaling is 
relevant for all driven thermal systems, including flow- 
ing liquids, mechanically deformed glasses, and stretched 
proteins. Quantitative analyses of these driven thermal 
systems, based on modifications of Eyring's theory to 
take the proper scalings of the barrier height into ac- 
count, will lead to an improved understanding and de- 
scription of these systems. 
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